Evaluation of the Permutational Structure of Quantum Gases at Finite Temperature 
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Proposed is an alternative method for permutational sampling in quantum gases using the path 
integral formulation of statistical mechanics. It is shown that in principle we are able to use two 
operators which enable us to construct a Markov chain through a graph of the irreducible represen- 
tation of the symmetric group. As an illustration of this method, a test calculation of four particles 
in a harmonic trap is performed. 



INTRODUCTION 

Irrespective of their experimental realisation, the the- 
oretical existence of condensed phases of noninteracting 
media is indeed remarkable, and the mathematical struc- 
tures with which they are imbued can lead to a variety 
of interesting observations. It is not difficult to show 
that a system of noninteracting indistinguishable parti- 
cles with bosonic symmetry will undergo a second order 
phase transition if the spatial dimension is greater than 
two. For d < 2 it has been rigorously been shown that 
a bosonic system will not undergo phase transition for 
T > [l|, Q • Theoretically however it was found that this 
is dependent on the form of the confining potential 3] con- 
firmed by evidence of Bose-Einstein condensates (EEC's) 
having recently been observed in quasi-2d traps ^j. 

The experimental observation of EEC's in systems 
with both positive and negative s-wave scattering 
length [^, 0,15, and the rapid expansion of novel trapping 
potentials in EEC's such as optical and magnetic lattices 
@i @i fipt and quasi-one and two dimensional harmonic 
traps (j, IllL Il2l |. has given a renewed impetus to the 



theoretical description of low temperature atomic gases. 

Of the theoretical approaches, path integral Monte 
Carlo(PIMC) techniques stand as one of the most useful 
non-perturbative methods, and is the only method able 
to produce exact properties of systems at finite tempera- 
ture, via sampling of the thermal density matrix. One of 
the ubiquitous problems in all approaches to evaluating 
the partition function for bosonic systems is a transpar- 
ent way of sampling the permutation space inherent in all 
problems involving indistinguishable particles. Previous 
approaches such as those in [13|, |l5| advocate either 
the sampling of some cyclic subset of Sn , the symmetric 
group of order n!, typically up to C4 or C5, or approxi- 
mate the interacting partition function structure to that 
of a noninteracting system, which is only strictly valid 
for systems which are weakly coupled. 

In this paper we wish to make some remarks regard- 
ing the permutational structure of the partition function 
for bosonic particles. We present an algorithm which al- 
lows us to construct a Markov chain through permutation 
space for an explicit probability measure which maps be- 
tween conjugacy classes of the symmetric group. It is 



partially inspired by the recent work regarding Markov 

m, [13, El, 111 and the con- 



processes on Young tableaux 
nection between irreducible representations of the sym- 
metric group and the analytic form of the bosonic parti- 
tion function for noninteracting systems. Unlike a previ- 
ous method ^15i] . which is restricted to weakly interact- 
ing systems, we propose that this method samples the 
ring configuration structure at any temperature without 
any a priori knowledge of the noninteracting counter- 
part. This enables us to quantify ring configuration prob- 
abilities over a larger temperature scale and evaluate the 
relative importance of certain structures to the partition 
function. 

The aims of this paper are (1) to give a rigorous formu- 
lation of the role of permutation cycles in the canonical 
ensemble partition function via a correspondence with ir- 
reducible representations of the symmetric group, and (2) 
to propose a practical scheme for computational imple- 
mentation within PIMC. As an illustration we perform 
calculations on noninteracting bosons in a 1-dimensional 
harmonic potential, the results of which are presented 
in sec. IV, a system which has a very well understood 
theoretical basis. 

This paper is organised as follows. In Sec. II we pro- 
vide the theoretical outline for the path integral formu- 
lation of statistical mechanics in the canonical ensemble 
and derive a form of the partition function based upon 
the conjugacy classes of the symmetric group of n ele- 
ments. In Sec. Ill we outline the basis for effective eval- 
uation of the discrete permutational space by considering 
random growth of partitions of the symmetric group us- 
ing one loop operators. We are able to stochastically 
sample the graph of partitions and show that this the- 
oretically will enable us to determine loop configuration 
probabilities as a function of temperature. 



PATH INTEGRAL STATISTICAL MECHANICS 

The basis of the path integral formulation of statisti- 
cal mechanics is the decomposition property of the many 
body density matrix p{R,R';f3), and the resulting clas- 
sical configuration integral for calculating the properties 
of quantum systems [isl . [20l | 
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piR,R';P)^ j dRip{R,Ri; 13/2) p{R^,R'- (3/2) (1) 



where 

= (i?|e-^«|i?') (2) 

and \R) = I ) denotes the coordinate repre- 

sentation of the fuU N particle Hilbert space. The Hamil- 
tonian operator TL in general contains q-number quan- 
tities representing particle kinetic and potential terms 
and Eq. (2) represents that of an interacting many 
body quantum problem for distinguishable (Boltzman- 
nian) particles. The partition function is given as the 
trace of the density matrix 

Z = Trip) 

dR{R\e~^^\R) 



Tr (e-^«) 



(3) 



the final line indicating that it is independent of the 
representation. The convolution property of the den- 
sity matrix is exact and gives us a representation of 
the density matrix at temperature T as the convolution 
of two density matrices at temperature 2T, effectively 
mapping low temperature problems to that of higher 
temperatures. Using the Trotter formula, exp {—f3H) = 
exp (— tH) , where t = /3/M, we may perform this op- 
eration an arbitrary number of times to obtain the ex- 
pression 



M 

n 

1=0 



dR,p{R,+i,R,-(3/M) 



(4) 



with the trace condition requiring that R^.j-^i — Ro- 
In the limit that M tends to infinity, t and we 
obtain an expression for the partition function which is 
functionally dependent on the classical action. The path 
integral expression is 
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For the discrete form of the partition function given 
in Eq.(4), the procedure of taking the trace of a convo- 
lution of M density matrices is to represent particles as 
a set of M classically interacting beads interacting via a 
temperature dependent harmonic force, forming a closed 
'necklace'. This is sometimes referred to as the polymer 
isomorphism It should be noted however that the 

inter-particle potential only acts between beads of cor- 
responding time-slices (M) [3, [23|. As the number of 
time slices tends to infinity particles are represented as 
continuous closed loops. For an extensive review of path 
integrals in this context, the reader is referred to [13[. 

The formal evaluation of Eql5] for noninteracting par- 
ticles, gives the partition function to be the product of 
the single particle partition function, i.e. 



Z 



(7) 



When considering ensembles of indistinguishable par- 
ticles one is required to permutate over particle labels. 
In the discretised form of the partition function we need 
only consider permutations over particles in the final time 
slice [isl, i.e. 



pb{R,R';P) - ]^ E p{R.<i>R'\P) 



Zb = 



n 



dR,p(R,+uR,;l3/M) 



dRMp{RQARM.l3/M). (8) 



For noninteracting systems of identical particles, the 
convolution property of the density matrix leads to a 
decomposition and the partition function becomes the 
product of density matrices at various temperatures. 
Consider a system of two particles with the Hamiltonian 
operator 



Ti. = Hi + H2 



(9) 



where 



such that [7^1,7^2] = 0. In this case the Boltzman- 
nian 2-body density matrix will be p{xiX2, x[x2', P) = 
p{xi, x'l, (3)p(x2, X2', (3). The partition function for the 
bosonic case will be 



IvR ^ lim / dR, 



(6) 



is the functional integration measure. The dynamical 
quantity of interest, the Euclidean action, is now ren- 
dered a c-number quantity on which the partition func- 
tionally is dependent. 



dxidx2p{xi,xi; I3)p{x2, X2; P) 



+ - J dxidx2p{xi,X2;l3)p{x2,xi;l3) (10) 

Using Eq.(l) for the second integral in Eq.(lO), the 
partition function can be written exactly as 
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Zsm = \ [z,[py + z,m] , (11) 

Zi denoting the single particle partition function. As 
is frequently noted, the partition function of statistical 
mechanics is related to the propagator of quantum theory 
via rotation of time to the imaginary axis [? ] , and so the 
effect of permuting particle labels in the two particle case 
is to wind the trajectory twice around imaginary time 

A general expression can be obtained for the N- 
body noninteracting partition function, involving the ir- 
reducible representations of the symmetric group. Con- 
sider the symmetric group, denoted by Sm which is the 
group of all permutations on n objects, with the order 
being \Sn\ = 'rL\. The number of conjugacy classes of 
Sn is equal to the number-theoretic partition function 
P(n) [21], with the property 




FIG. 1: Representation of the Conjugacy Classes of Si 



P(n) 
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(12) 



Here 0' is the i*'' conjugacy class. The conjugacy 
classes can be graphically represented via the standard 
Young tableaux, however for this case it is instructive to 
consider diagrams as shown in Fig.(l). This originates 
from the representation of particles in path integral sta- 
tistical mechanics as kinetic strings. The action of el- 
ements of Sn on these strings is to attach each string 
starting position (r — 0) with all possible final positions 
and integrate, giving the partition function. It should be 
noted that this is exactly true for the interacting as well 
as the noninteracting case. 

Consider 5*4 which contains a total of 5 partitions and 
24 elements, which is shown in Fig 1. From this is can 
be seen that elements of a particular conjugacy class are 
topologically equivalent. We can identify Young tableaux 
and partitions of n with the weakly decreasing sequence 



0^ = {Al,A^,...,Al}, 



each element having the property that A^^ > A2 > 
Ap, or as the strongly decreasing sequence 



(13) 



> 



0' = {{x\)'-\{x^) 



where X\ > \i > ... > Xt. 



(a;)"'-} (14) 



Here kij will be denoted as 
the multiplicity of a particular loop size. The dimension 
of each conjugacy class is equal to the number of Young 
tableaux with shape (jj^ which is given by 



n! 



nu ((^ 



k I 



(15) 



For the purpose of this paper we need to identify two 
numbers |(/)*| — p, the number of elements in a partition, 
which is equal to the number of rows on the equivalent 
Young tableaux and = 7, the number of distinct 
elements. For example if we denote $ as the set of all 
Young tableaux with n boxes, for the case of S4 we have 
the conjugacy classes represented by the partitions 



$ = {{1,1,1,1},{2,1,1},{2,2},{3,1},{4}} 
= {{in, {2, in, {2^, {3,1}, {4}} 

= {0^0^0^0^0n, (i6) 

with the dimension of each conjugacy class given as 



dim ($) = { 



24 



24 24 24 24 



} 



14.4!' 2.12.2!' 22.2' 3.1' 4 
= {1,6,3,8,6}. (17) 

With this notation, we are able to write down a general 
expression for the partition function of a noninteracting 
bosonic system of n particles, which is 



Zf 



. Pi") / i0i \ 

dim(0On^i[^}/5]'''M 

i=l \ j=l J 



(18) 



where the superscript ki^ denotes raising the power of 
the partition function to the multiplicity of this particular 
loop with cycle size A*. To obtain an understanding of 
Eq.(18), let us again consider the example of S4 for which 
P(4) = 5. Written out fully, the partition function takes 
the form 
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FIG. 2: (Color online). Loop structure of the 4 particle par- 
tition function in a 1-d harmonic potential. The x-axis is the 
temperature and the y-axis the weight of the i*'^ conjugacy 
class. 



I0'l 

24xZb ^1[Z, [X](3] + 6 n ^1 [A?/3] + 
i=i i=i 



611^1 m"'^' 



6Zi[2/3]Zi[/3]2 + 3Zi[2/3]^ 



+8Zi[3/3]Zi[/3] +6Zi[4/3] 



(19) 



This is the exact form of the bosonic partition func- 
tion for 4 noninteracting particles in an arbitrary exter- 
nal potential. The probability of randomly choosing a 
particular partition 0* g $ as a function of temperature 
is given by 



, 10.1 
f^i^^e^) = — dim {cj,^) II Z,[Xlpf-^ 



(20) 



An analysis of this structure can show us the relative 
weight of partition structures at various temperatures. A 
graph of the contributions of the partitions of 4 to the 
partition function for noninteracting particles in a 1-d 
harmonic trap is shown in Fig 2. 

In this figure, the contributions from various elements 
of the permutation group are clearly evident, showing 



that at high temperatures the identity dominates, leading 
to the usual classical statistics 13|, |20[ . Also in the T ^ 
limit, the partition function reduces to 



ZbW- 



P{n) 

2 = 1 



(21) 



which is exact for all noninteracting bosonic systems. 
For systems in which the particles interact, the evaluation 
of partition function elements will generally not be de- 
composable as a product of single particle functions , 
but the loop structure based upon the conjugacy classes 
of Sn will still be present. However the range of prob- 
lems for which Eq.(8) is exactly solvable for interacting 
systems is severely limited. 



PIMC AND PERMUTATIONAL SAMPLING 

The PIMC technique at present is the only method 
which renders Eq.(8) amenable to numerical solution. 
However the method by which the permutational struc- 
ture is accounted for either is based upon including only a 
subset of cyclic exchanges, usually up to C4 or C5, or the 
approximation of the exact permutational structure with 
that of a noninteracting system. The first method works 
very well in predicting properties such as the superfluid 
transition temperature of helium [isl [2^ . but is equiva- 
lent to the treatment of parastatistics in bosonic systems 
first outlined by Green [22 



24| 



and does not include the 
full permutational structure of the system. The second 
method works only for weakly interacting systems around 
the superfluid transition point as implemented in ^15.] via 
knowledge of the non-interacting partition function coun- 
terpart and is unable to give describe systems where the 
interactions are strong. The strong interaction regime is 
important especially for the case of repulsive interactions. 
It is specifically this regime in which the PIMC method 
is superior as perturbation theory becomes divergent. 

In recent years mathematical methods have been de- 
veloped regarding the random growth of partitions of the 
symmetric group, the probability measure known as the 
Plancherel measure and random matrices 
One would like to ask if this abstract mathematical 
method could help in the very physical application of 
predicting the permutational structure of EEC's and su- 
perfluids at arbitrary temperatures. If we denote $ as 
the set of all Young tableaux with n boxes, the probabil- 
ity of randomly choosing a particular partition 0* e $ as 
a function of temperature is given by 



e$) = — dim((/)^)Z[4;/3] (22) 

with Zb equal to the total interacting bosonic partition 
function. This is presumably the probability used in [l5| 
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to construct partitions, which was then apphed to the 
case of a weakly interacting hose gas, however the way 
they did this was not exphcitly stated. 

In the simplest implementation of a stochastic sam- 
pling method, one would randomly select a partition and 
select this with the probability given above. However, as 
derived by Ramanajan and Hardy, for large n the number 
of partitions is given by [2^ 



P{n) 



1 



(23) 



4\/3n 

which becomes prohibitively inefhecient for large n, as 
the rejection rate for transitions between uncorrelated 
partitions would be large. 

We propose a set of operators which enable us to con- 
struct a Markov chain through the conjugacy classes of 
Sn- a similar formalism can be found in Borodin and 



Olshanski [17|, ll8| in constructing partitions, identifying 
an operation which maps between partitions. Again con- 
sider a partition as the weakly decreasing sequence 



{A^, A2, A*, A^_;^, A* } 



(24) 



but with the restriction that at least one element, 
namely Ap, has the value of 1. Consider two operators aj 
and cij , such that 

a,q^^ = {ai,...,(a} + i),...,a;_j 



"1' 

{2,13} 
{2^,1} {3,12} 

{3,2} {4,1} 



{5} 



FIG. 3: Neighboring partitions of under the action of a 
and a 



One can easily see that these transition probabilities 
satisfy the criterion that the sum of probabilities over 
neighbours is equal to one. In the construction of a 
Markov chain through permutation space, one could sup- 
pose that these are the correct transition probabilities to 
be used. Ergodicity, in the sense that the correct weight 
of each partition will be reproduced, is not assured via 
these relations. It is not however difficult to prove that 
all partitions are accessible under this scheme. If we con- 
sider the partition denoting the identity of Sn as {l}Ar, 
then we may construct any other partition using a re- 
peated application of aj , that is 



{Ai, A2, Ap} = n^^'t^j'^- 



(28) 



5,0^ = {ai,...,(a;-i),...,a;,a;+j 



(25) 



such that 



\aj(l),\ 



1. 



(26) 



These operators create and destroy Ci subgroups of the 
conjugacy classes, in the process creating a new partition 
which is also a conjugacy class of the relevant permuta- 
tion group. The action of these operators is to effectively 
move up and down a graph with P{n) vertices of the 
partitions of Sn , leaving the sum of the elements of par- 
titions invariant (cf. fig. 3). 

If we denote partitions which are related via a single 
application of aj or aj as neighouring partitions, we can 
construct a probability measure over neighbouring par- 
titions to give a formula for the transition probability 



under aj and aj as 



dim {aj4>i) 
J2j dim (aj (/)j) 

dim (qj (/).,) 
E, dim(aj 0j) 



(27) 



Since we are able to perform the reverse of this oper- 
ation, that is reach the identity from any partition, via 
the use of aj, all partitions are connected. Once a parti- 
tion is chosen then it can be accepted or rejected via the 
Metropolis scheme [26l|. One hinderance to the successful 
application of these transition probabilities is that there 
does not exist any formula for the sum over neighbours 
which would be relatively easy to implement on a com- 
puter within a Monte Carlo code, especially for systems 
containing a large number of particles. 

An alternative option which is relatively easy to imple- 
ment would be to give all neighbours the same weight. 
Then one may choose a particular element of a partition 
to act upon with the probabilities 



{(l)i;aj(j),) = 



10.1-1 



10. 



(29) 
(30) 
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The action of aj on a partition element is to give 
a(A) (A — 1, 1) and as such the normalisation is the 
number of elements in a partition which are not equal to 
one. The action of aj on a partition element is to give 
a(A, 1) — > (A-l-1) such that at least one element must have 
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value one, and as such the normahsation is the number 
of elements in a partition less one. 

It was essentially this use of aj exclusively by Bonin- 
segni in to construct permutation cycles for the su- 
perfluid phase of ^He, containing 64 atoms. However in 
this work permutation cycles were constructed at each 
sampling move from the identity. Again, from Eq.(28) 
we can see that each permutation cycle will be accessible 
via this scheme, however it is not guaranteed that the 
correct weights will be reproduced. 



COMPUTATION AND DISCUSSION OF 
RESULTS 

As a test of the proposed scheme with the analytic 
results given in Fig. 2, we performed a path integral 
Monte Carlo calculation on 4 bosons in a 1-d harmonic 
trap. We used essentially the method as outlined in[l^, 
in the primitive approximation 27j . 

The Trotter number was chosen such that the single 
particle energy was comparable to the analytic expression 
{E) = icoth(i/3), however small enough that chosen 
permutations were accepted. 

We implemented two schemes for permutation sam- 
pling. In the first, the results of which are summarised 
in Fig. (4), particle labels were randomly shuffled at 
each move and were accepted or rejected according to 
the Metropolis scheme. This was used primarily as a test 
for the validation of the code. In the second scheme we 
constructed a random walk through permutation space 
based upon the transition probabilities given in Eqs. (29- 
30) between partitions. At each move a partition list is 
created and an up or down move {aj or cLj) proposed 
with equal probability. In the case that a partition is 
at the end of a path (cf. Fig. (3)), then the algorithm 
forces a move to the next connected vertex. In the case 
of Qj , once a partition element is chosen, then the nearest 
single loop structure is found such that the acceptance 
probability in the following Metropolis sampling move is 
maximised. The results of this method are shown in Fig. 
5. 

As can be seen, randomly shuffling particle labels re- 
produces the analytic decomposition of the partition 
function very well. We would however expect this 
method to become diminishingly inefficient as the par- 
ticle number increased, and as the confining potential 
becomes significantly weaker, where the thermal wave- 
length of the particles will no longer be comparable to 
the mean particle separation, and thus would not be a vi- 
able option for a general PIMC code. Further the results 
of our method show that even though every partition 
is accessible via this algorithm, this does not guarantee 
that the method will reproduce the analytical result of 
Eq.(20). 

An interesting feature of the results of the two meth- 




FIG. 4: (Color online). Partition structure of four bosons in 
a 1-d harmonic trap from PIMC calculations, with partitions 
created via a random shuffling of particle labels. The solid 
line is the exact answer as given by Eq.(20) 




T/kB 



FIG. 5: (Color online). Partition structure of four bosons in 
a 1-d harmonic trap with partitions created the new method. 
The solid line is the exact answer as given by Eq.(20) 
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+ Random Shuffle 
^ Guided sampling 



FIG. 6: (Color online). Energy as a function of temperature 
for 4-particles bosonic particles using two different methods 
for sampling permutation space. The statistical errors in these 
calculations are beyond the resolution of this graph. 



However the overhead required for the construction of 
partitions is minimal, with only knowledge of the cur- 
rent permutational structure and the proposed structure 
required, and thus memory requirements kept to a mini- 
mum. Further it is not necessary that the entire graph of 
partitions be sampled as many configurations will have 
only minor contributions, these depending on the tem- 
perature. After an initial period of equilibration, the 
process will sample the local neighbourhood of the high- 
est weight permutations, those which will contribute the 
greatest to observables. Although the intrinsic n\ scal- 
ing involved in the sampling of permutation space has 
largely been overcome, they have been replaced by al- 
gorithms scaling exponentially 13, 14, 3], which is still 
problematic for large n. A recent proposal for a size inde- 
pendent algorithm has been reported by Boninsegni et. 
al. in [2^ which also looks promising for investigations 
into large systems and calculating off-diagonal elements 
of the density matrix. 



CONCLUSION 



ods of permutational sampling investigated is that, whilst 
they did not coincide with regard to the predicted per- 
mutation cycle structure, the resultant energy of the two 
methods was nearly identical (cf. Fig. (6)). 

Although the new method proposed was unable to re- 
produce the exact analytic form of Eq. (20) for this rela- 
tively basic model, we would expect that as the number 
of particles is increased that this would become more ac- 
curate. As the number of paths between vertices of the 
partition graph is increased, a walk on this graph will not 
be required to pass through a disproportionate number 
of low weight vertices. Our results show that the weight 
of the larger cycle lengths is underestimated, which is 
reflected by the fact that they lie on the end points of 
the graph of 5*4. The larger cycle lengths contribute to 
quantities such as the mean winding number which de- 
rived from the mean squared winding number (2^. This 
method does have the advantage over other methods in 
that the only approximation is the form of the transition 
probabilities. After a period of equilibration, where the 
cycle structure is constructed from the identity, we will 
be able to sample the local permutational structure by 
using aj and Oj without any a priori knowledge of cyclic 
structure of the partition function. 

As mentioned earlier, the number of conjugacy classes 
of Sn has the asymptotic form 



P{n) 



4\/3n 



(31) 



To summarize, we suggest the use of one loop operators 
acting upon partitions for the construction of a Markov 
chain through permutation space. We have derived a 
form of the partition function for bosonic systems in the 
canonical ensemble based upon the irreducible represen- 
tations of the symmetric group of n objects. It was shown 
that the partition function can be decomposed into a sum 
over the conjugacy classes of Sn which defines a proba- 
bility measure over Young tableaux which is a function 
of temperature. 

There have been various proposals for a more trans- 
parent method of sampling permutation space in the re- 
cent literature both for bosonic and fermionic systems 
13, 3 28, 2^, 13], which appear to be successful in pre- 



which would provide the associated computational 
sampling frequency as a function of particle number. 



dieting properties of interacting quantum systems. How- 
ever there still does not exist an unambiguous method for 
the sampling of permutation space at the same level of 
rigor as for calculating expectation values of continuous 
probability measures. One may expect that a multilevel 
Metropolis scheme could be developed whereby the par- 
tition structure at finite temperature is reproduced ex- 
actly, giving a more accurate sampling of the partition 
function. Any scheme that purports to full sampling of 
the permutational structure of atomic gases within path 
integral Monte Carlo should be able to reproduce the an- 
alytic form for the non interacting case. A more rigorous 
formulation for the construction of discrete probability 
measures amenable to Monte Carlo methods may lead to 
some very interesting new observations. 

An interesting mathematical question also arises which 
we wish not to explicitly address here, but may be a use- 
ful connection between the theory of indistinguishable 
many-body systems and a large body of work in probabil- 
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ity theory and the random growth of Young tableaux. Is 
one able to construct a Markov chain via the use of aj and 
CLj in the canonical ensemble that in the infinite time limit 
will give probabilities distributed according to equation 
20? Further, what is the limit shape of Young tableaux 
in the N oo limit as a function of temperature [l^, [3]? 
As r ^ oo we expect this to be x9 (x + where 9 is 
the Heaviside step function, as the identity will become 
dominate in this limit. In the T limit, where the 
partition function reduces to Eq.(21), the limit shape in 
the grand canonical ensemble has been established and as 
such one would expect a continuous transformation be- 
tween these limit shapes, the mapping being dependant 
on temperature. This will be related to finding eigenval- 
ues of random matrices as a function of temperature. We 
hope to investigate this further in a future paper. 
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